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Abstract: Neutron computed tomography (NCT) is widely used as a noninvasive measurement 
technique in nuclear engineering, thermal hydraulics, and cultural heritage. The neutron source 
intensity of NCT is usually low and the scan time is long, resulting in a projection image containing 
severe noise. To reduce the scanning time and increase the image reconstruction quality, an effective 
reconstruction algorithm must be selected. In CT image reconstruction, the reconstruction algorithms 
used can be divided into three categories: analytical algorithms, iterative algorithms, and deep learning. 
Because the analytical algorithm requires complete projection data, it is not suitable for reconstruction 
in harsh environments, such as strong radiation, high temperature, and high pressure. Deep learning 
requires large amounts of data and complex models, which cannot be easily deployed, as well as has a 
high computational complexity and poor interpretability. Therefore, this paper proposes the OS-SART- 
PDTV iterative algorithm, which uses the ordered subset simultaneous algebraic reconstruction 
technique (OS-SART) algorithm to reconstruct the image and the first-order primal-dual algorithm to 
solve the total variation (PDTV), for sparse-view NCT three-dimensional reconstruction. The novel 
algorithm was compared with other algorithms (FBP, OS-SART-TV, OS-SART-AwTV, and OS-SART- 
FGPTV) by simulating the experimental data and actual neutron projection experiments. The 
reconstruction results demonstrate that the proposed algorithm outperforms the FBP, OS-SART-TV, 
OS-SART-AwTV, and OS-SART-FGPTV algorithms in terms of preserving edge structure, denoising, 
and suppressing artifacts. 
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I INTRODUCTION 


Neutron tomography is a noninvasive measurement method that differs from X-ray imaging, and 
muon tomography has been applied in various fields with a high accuracy and reliability [1-2]. The 
strong transmission of neutrons in metallic elements makes neutron tomography extremely useful for 
examining hydrogen in heavy elements under harsh experimental conditions such as strong radiation, 
high-temperature, and high-pressure two-phase flow systems. In addition, neutrons are used in boron 
neutron therapy owing to their ability to detect and reconstruct the three-dimensional (3D) distribution 
of boron concentration [3-4]. Neutron tomography systems consist of two parts: a neutron photographic 
system and a computed tomography (CT) system capable of measuring and displaying the 3D void rate 


in the boiling flow in a heated rod beam. They have also been applied to simulate the cores of advanced 
nuclear [5-8]. Asano and Takenaka [9] used neutron tomography to determine the void-rate distribution 
in two-phase air and water flows. Kureta et al. (2008) visualized the toroidal streaming in a heated fuel 
rod bundle using 3D neutron tomography [10]. Tremsin et al. (2013) used neutron imaging to assess the 
attachment, deposition, and structural integrity of fuel core blocks [11]. Andersson et al. (2015) 
developed a portable fast NCT to obtain the void rate distribution in fuel rod bundles in a boiling water 
reactor [12]. Although portable NCT systems can perform measurements at any altitude, their neutron 
generators exhibit extremely low neutron yields. Therefore, promoting the rapid development of 
neutron tomography for improving the efficiency and safety of nuclear fuels is an important research 
direction. 

The most crucial part of a neutron tomography system is the 3D reconstruction algorithm, which 
is directly related to the image quality and measurement reliability. Because of the long scanning time 
of neutron CT (ranging from 30 min to several hours) and the complexity of practical application 
environments, it is difficult to obtain complete projection data. Consequently, reconstruction algorithms 
that can efficiently utilize sparse-view projection data are required to obtain accurate measurement 
results. Furthermore, the reconstruction of sparse-view projection data is helpful for investigations in 
harsh environments such as strong radiation. Currently, algorithms are used for the reconstruction of 
sparse views and highly under-sampled projection data. Therefore, this study was focused on a 3D 
reconstruction algorithm for sparse-view NCT. 

Sparse-view reconstruction is a highly complex mathematical problem and is a significant 
research direction in the image reconstruction field. Because prior information can be embedded in the 
iterative reconstruction algorithm, it is considered the sparse-view reconstruction algorithm most likely 
to replace the FBP algorithm. With the development of deep learning, it has been applied to the field of 
tomographic reconstruction. Many researchers have worked on sparse-view tomography and achieved 
significant results. As proposed in [13], the DEAR model, which adds a priori information and data in 
the image domain to the compressed sensing-based variation model, can not only eliminate the artifacts 
of the reconstructed image but also improve the quality of the reconstructed image. A multiple 
adversarial learning angiography image reconstruction framework has been proposed in the literature 
[14], which solves the challenge of low-intensity aortic reconstruction by introducing dual-correlation 
constrained adversarial learning; its application in clinical data shows its feasibility and effectiveness. 
However, because of practical limitations, this study focused on neutron CT 3D reconstruction, and it 
was not possible to obtain a large number of neutron projection images for supervised training. 
Therefore, this study focuses on an iterative algorithm for neutron CT 3D reconstruction. Currently, 
mainstream iterative reconstruction algorithms include ART [15], SIRT [16], and SART [17]. Although 
all these iterative algorithms can reconstruct high-quality images, when the projection data are 
extremely sparse and no a priori information is introduced, the reconstructed images show significant 
artifacts. Donoho [18] proposed the compressive sensing (CS) theory, which has promoted the rapid 
development of sparse-view image reconstruction. According to the mathematical implications of the 
CS theory, an image may be reconstructed exactly if it is sparse or can be sparsely represented by a 
spatial transformation. The evolution of the CS theory has facilitated the application of total variational 
algorithms in sparse-view image reconstruction [19]. Rudin proposed the ROF denoising model, which 
reduces noise while preserving the image edge and detailed structural information using the total 
variance as a constraint for image denoising [20]. Several TV-based techniques and variations to 
improve the performance of the algorithm have been put forward, including directional TV [21], 
weighted TV [22], edge-preserving TV [23], and weight TV[24]. The total variation model is a typical 
regularization model. Total variance minimization is particularly important in CT imaging because of 
its ability to resolve acute discontinuities. This is critical for many imaging problems, because the 
sample edges contain critical structural information about the object. However, because of the non- 
smooth character of the total variance, it is difficult to achieve total variance regularization 
minimization. Furthermore, one drawback of these models is that they assume that the image is 
piecewise constant, which destroys significant information in the image. The total-variational image- 
denoising problem is a class of convex optimization problems. To solve convex optimization problems, 


the second-order method exhibits good convergence and requires only a small number of iterations. 
However, each iteration of the second-order method is extremely complicated, making it difficult to 
apply to large-scale problems. In contrast, the first-order method involves only function values and 
gradient information, and the process of each iteration is relatively simple. In the case of gradient 
magnitude sparsity, the total variance of the intermediate images is minimized, resulting in high-quality 
reconstructed images. Therefore, this study employs a first-order method to address the total-variance 
image noise reduction problem, and proposes an algorithm applicable to the 3D reconstruction of NCT 
in sparse views. 

The goal of this study was to find an efficient 3D reconstruction algorithm for sparse-view NCT. 
To address artifacts and noise in sparse-view projection-reconstructed images, we propose the ordered 
subset simultaneous algebraic reconstruction technique (OS-SART-PDTV) algorithm, which uses the 
ordered subset simultaneous algebraic reconstruction technique (OS-SART) to reconstruct the image 
and the first-order primal-dual algorithm to solve the total variation (PDTV). This novel algorithm 
achieves a blurred version of a piecewise constant object through phenomenological modeling, and 
reconstructs high-quality 3D images by minimizing the total variance of the intermediate images 
through gradient-amplitude sparsity, which allows the reconstructed image interiors to change quickly 
and smoothly. Next, we present the proposed method, including the proposed algorithm and NCT 
fundamentals, in Section 2. In Section 3, the simulation experiments and a real neutron projection 
experiment are discussed. Section 4 discusses the research process and presents conclusions. Finally, a 
short conclusion is presented in Section 5. 


I PRINCIPLES AND ALGORITHMS OF NCT 


2.1 NCT system 


This section describes the fundamentals of CT image reconstruction. In general, CT image 
reconstruction generates a 3D volume from a 2D collection of projected images using mathematical 
algorithms. This requires the development of a measurement model using mathematical methods to 
relate the measured data to the desired physical properties. Image reconstruction consists of two 
processes: forward and back projection. Forward projection enables the calculation of model 
measurements corresponding to the actual measurements that match the physical properties. Model 
measurements corresponding to the actual measurements were calculated using forward projection. 
Forward projection refers to the mapping of an objective onto the projection domain. Back projection is 
the opposite of forward projection, which determines the physical properties of a measurement. 
Therefore, these two operations largely determine the image reconstruction precision. The change in 
forward projection depends on the radiation beam geometry. Based on the selected system, source of 
radiation, and available data acquisition systems, the appropriate geometry of the beam was selected so 
that reliable measurements could be obtained. The NCT system primarily includes a collimator and 
shield, rotating sample stage, and neutron imaging system. The neutron imaging system included a 2D 
neutron image detector, conversion screen, and 3D reconstruction software. The algorithm used for 3D 
reconstruction directly determines the imaging quality and measurement reliability of the system, and 
is the core of the system. A schematic of the NCT scanning system is shown in Figure 1. 
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Fig. 1 Schematic of the NCT system 
In NCT systems, the mass of the neutron source and noise generated by the electronics can cause 
the projection data to contain noise. Therefore, Poisson noise and Gaussian noise are added to the 
numerical projection data reconstruction to validate the denoising ability of the algorithm while testing 
its reliability [25]. Therefore, these two noise components must be considered during NCT. The 
equation is as follows: 


1 . = Poisson(I,) + Gaussian(m,, + 82) ) 
where I denotes noisy transmission data, I , represents the mean photon number, M ; o 
indicates the average value of the electronic noise, and 6 2 is the variance of the electronic 
noise. System calibration typically has a value fromm ; o to zero. The variance 6 2 was 


estimated via dark current measurements [26]. 
2.2 Principle of CT 
2.2.1 CT projection acquisition 


In NCT, the geometry of the neutron beam is primarily parallel. Among the mathematical concepts, 
the parallel geometric beam is the simplest. An object can be represented by a binary function R(x, y), 


which represents the intensity function on the x-y plane. The projection represents the decay integral 


of the neutron beam as it passes through the object from the neutron source, which is considered to be 
the sum of the line integral and ray. The expression of the projection is given by the equation 


h(s, B) =| Rœ, y)ds Q) 


where L represents a neutron beam passing straight through an object R(x, y), S represents the 


detector coordinates in the projection domain, and 2 denotes the counterclockwise rotation angle of 


the x-axis in the projection domain. According to the Delta functions in Eq. (2) can be rewritten as 
follows: 


h(s,B)=[_ | R(x, y)d(xcos B+ ysin B—s)dxdy 8) 
Eq. (3) represents the Radon transform of the object R(x, y) , which transforms the space data into 


the frequency domain. The projection consisted of a group of integrals of the neutron attenuation 
coefficients of the sample measured at the detector for every given angle. 


2.2.2 Central slice theorem 


The Central slice theorem establishes a connection between the 2D projection of an object and its 
1D projection in the Fourier domain. The 1D Fourier transform function K(o) of the projection 


function A(s, 8) of the density function R(x, y) in a certain direction is the value of the 2D Fourier 
transform of the density function R(x, y) on the x-y plane along the same direction on a straight 


line past the origin. Eq. (4) represents the defined 1D Fourier variation of the objective function in 
polar coordinates. 


K(o, B) = K B R(x, yje ee dydy (4) 


Eq. (4) represents the 2D Fourier transform of the spatial frequency ox =ocosf,oy=osin# in the 


Fourier space. Therefore, the following equation was derived: 


2m po f E TEE 
R(x, y) =Í f, O(o,,0, Je? EE y dod B (5) 


Several problems arise when the central slice theorem is used for the Fourier transform 
reconstruction. First, the theorem produces data in the Fourier space that are inconsistent with the 
Cartesian space. The data were interpolated into Cartesian coordinates. However, interpolation in the 
Fourier space before the anti-Fourier transform can have a significant influence on the reconstruction. 
The second problem is the difficulty in performing targeted reconstruction and identifying fine 
structures within small areas. 


2.3 Algorithms 
2.3.1 TV algorithm 


In the NCT, the measurement time can be drastically reduced by reducing the number of projection 
views; however, this leads to sparse sampling. Based on the CS theory, a signal can be reconstructed 
exactly if it is sparse or can be represented sparsely. Under ideal conditions, NCT reconstruction can be 
converted to solve the problem of a linear system; thus, the discrete model of CS theory may be 
represented. 


Lg =b (6) 
where L= {1 


;, f stands for system matrix b denotes the projection data obtained at each angle and g 


represents the image to be reconstructed. Mathematically, solving Eq. (6) usually converts it into a least 
square problem: 


arg min [Ze -bl (7) 


However, in sparse-view image reconstruction, the incompleteness of the projection data leads to 
reconstructed images containing noise and artifacts. Therefore, various regularization models have 
been proposed to obtain approximate solutions, which are defined by the following equations: 


: 2 
arg min |Lg— bl, +AZ(g) (8) 
. 2 
where arg min|Lg -b| 5 Tepresents the data fidelity term, AZ(g) denotes the regular term, and 
g 


A is the regularization parameter used to balance the regular and data terms. The equation that 
defines TV as regularization is given in Eq. (9) 


lel =DylV.al +y, e] +e ©) 
where g denotes the 3D image, V,g represents the gradient along i direction, and [e| denotes the 


complex modulus. 
2.3.2 First-order primal-dual algorithm 


Owing to the segmental smoothness of the regularization term of the variational model, the total 


variational model can effectively preserve image edges during image recovery and achieve improved 
recovery results. Therefore, the total variational model, particularly the total variational regularization 
model [27] proposed by Rudin, Osher, and Fatemi, has gained popularity. According to the ROF model, 
the discrete total variational model of the 3D images can be formulated as follows [28]: 


min +|Hg -uÈ +aCTV(g) (10) 
gs 2 2 


where |l*||2 stands for the Euclidean parametrization, a > 0 is the regularization argument, and 
CTV (g) is the regular term of the discrete total variational model. To obtain a stable solution to the 


problem of minimization, V , which denotes the difference algorithm on the space R3” a is defined 
as shown in Eq. (11) below: 


(Volos, = (VE) eng ’ (V2) ic 2 (V2) eng) (11) 


where 


Ee+l,u,q g Sewg e< N 


(VE) eug = | (e=0,1,2,..., N) 


0, e=N 
v Ee urq z Seu u < M 
(V2) .u¢ = (u =0,1,2,...,M) 
rae 0, u=M 
w Seutlq > Seng q < Z 
(V2) eu = (q =0,1,2,...,Z) 
na ‘ a i 


Zeu g denotes the ((N*M)*q+e+N*u) -th element of g , whose position in the image is 


(e,u,q); (VE) bug ; (V2) ug , and (VE) ong are the different operators in the x , y , and 


Zz directions at a pixel ((N * M )* q +e+ N *u) . The expression of the regular terms of the total 


variational model is shown below. 


crv(g)=|vel,= © SEOD 
O<e,u,gs 


In this study, the above model was improved using the finite difference matrix R instead of the 


Í 


Oaa, HOD 


(12) 


difference operator V. 
CTV (2) =||Rell, (13) 


Thus, the expression of the model investigated in the present study can be obtained as follows: 
. 1 2 
minz |Hg -u|, +a [Re], (14) 


Although the total variational model efficiently preserves image edges, it is difficult to obtain the 
minimum value of the total variational model by solving Eq. (14) directly because of the non-smooth 
nature of the total variational model [29]. Because the first-order primal-dual algorithm is an efficient 
method for solving non-smooth convex optimizations of images, we employed it to solve Eq. (14). 

Definition 1: Assume that the binary function J(y,v)(y¢N",veN") is convex for y, and concave 


for v . If, for any (ye N",v e N”), there exists an (vv) such that Eq. (15) holds, (vv) is the 
saddle point of the function J(y,v). 
TOV) SI VSI) (15) 


Chambolle and Pock suggested a first-order primal-dual algorithm to solve the following model based 
on the properties of the saddle points [30]: 


min max £(7,v) = 0(7)+(K z,v)—- HV) (16) 


CN" 


where O and Q are lower semi-continuous convex functions, and K denotes the line operator. The 
first-order primal-dual iterations were as follows [31,32]. 


+ i i 
x =argmin WHK xv) -z-z k 
g = VY 4 (yO _ y) (18) 
yar) = arg min KZ a) y HV) 1 y-y® 2 (19) 
vcN” 2t ? 


where the parameter s,t>0 is the iteration steps for the primal and dual variables, separately; 0 
denotes the combinational parameter, which is used to ensure the algorithm convergence. 

Definition 2: Assuming that function o(g):N” —>N , we call the formula in Eq. (20) as a 
conjugate function of the function o(g) [33]. 


o'(Ẹ)= sup {57g -o(2)) (20) 


where & represents the dual vector of g, N denotes the real number space, N” represents the k- 


dimensional real space, and T denotes the transpose. 


T 
max 


Izel jas Re 
| can be expressed as FIs . Then, Eq. (14) can be 


According to Definition 2, the 


converted into a minimal—extreme ape [34] as follows: 


min max [Hg ulf, +aé"Rg (21) 
g |él.<t 2 
The iterative formula in Eq. (21) can be obtained by using a first-order primal-dual algorithm. 
1 2 1 2 
(d+) _ + _ TpT dd), + Jeo 
gi = arg max > Hg -u| + ag” R E” + 5 le g |, (22) 
gcN S 
z (d+) _ gt) 4 Og? ~ gi) (23) 
2 
d+1 d+ (d) 
EO =argmaxa E" Rg) é- l (24) 


The minimization problem for the subproblem g can be converted into solving the solution of the 


linear Eq. (25). H isa 3 x 3 a chunking matrix. 


1 1 
(H’H+-D)'g@ =(H’ut+—g© -aR E”) (25) 
S S 
where H T is the transpose of H i R 7 is the transpose of R 
1 
g“™ =(H"H +- ry (H'u+=g®° -aR"E®) (26) 
mg“ 


For subproblem € , the maximum Metn of € isfirst o into a minimization problem. 


Ea) = = arg max gé” R g g (d+) ses, =a _ Lje- (Eo — taR g Gi) a 


le. <1 


eH a 2 


Then, the equation Eq. (27) is the &“ =P,(é -taR g “*™) , which denotes the projection of 


E -taR g “* onto set A 


A= G eR” xR” :él < 1) (28) 


The defining equations for the projection operator of any vector q onset A are given by Eqs. (29) 


2 
2 > Vq (29) 


P,(q) =argmin|é—4| 


q 


Then, there is P,(q) = 
max(l, 


al) 


Table1 First-order primal-dual algorithm 


Set $ >0,t >0,8 €[0,1],input,u,a,R 


Initialize € ,¢ 


Repeat k<—k+1 


g“™ =H sty (H'u g -aR é”) 
S S 


gz (d+1) = gon +g - g”) 


Ee =P (é -taR z (a4) 


(i+1) (i) 
k” -e"l 


Calculate the deviation value RD RD = | Ey 
g 


2 


Until RD<10~, ork reaches the maximum number of iterations 


Output recovery image g¢“*” 


where s is the iteration step of the original variable, ¢ is the iteration step of the dyadic variable, and k is 
the number of iterations. 


2.3.3 Adaptive weighted total variation algorithm 


The conventional TV term is based on the assumption that the reconstructed image is piecewise 
constantly distributed, and this approach can cause excessive smoothing of the reconstructed image 
edges. To alleviate the problem of over-smoothing of traditional TV edges, many researchers have 
studied weighted adaptive total variation [35]; therefore, in this study, an adaptive weighted total 
variation (AwTV) minimization image model was used for comparison with the proposed algorithm. 

min |||, Subject to |p— Au|<e (30) 


u20 


where 4Awry denotes the adaptive weight of the total variance of the reconstructed image. The 


|a 


ay B defined as follows: 


2 2 2 
AwTV = > J, (Hy 52 fe + W, Hy yx =i yaz) T wW, (Uyy Yh yz) (31) 


XPZ 


la 


Hy y,2 = Hx-,y,z Hy y,2 a Hy y-1,z 


w, = exp| —( Y |,w, =exp]| —( y 
(32) 
y 


Hy yz 7 Mey cA 
w, = exp| (E2 ar 
T 


where w, denotes the weight factor, and T is a scaling factor in the weights used to control the 


diffusion intensity of each iteration. 


2.3.3 Fast Gradient Projection Algorithm 


Beck and Teboulle proposed a fast gradient projection algorithm [36], which is derived as follows. 


The model expression to solve the TV-based denoising problem is given by Eq. (33). 


min |x -b I +2ATV (x ) (33) 
EC F 


x 


TV is a nonsmooth regularization function that can be an isotropic T V z OF anisotropic 


TV ı -and C belongs to a closed convex subset of E =R ™ *" | The 


nonsmoothness of the TV function leads to significant difficulties in solving Eq. (33) .To address this 
challenge, Chambolle proposed a gradient-based algorithm for solving pairwise problems using a 
pairwise approach under unconstrained conditions. Thus, a constrained problem pair was constructed 


according to the proposed method. The following symbols are required. 


1 


© P denotes the set of matrix-pairs (p a) 2) which meet: 


(p : 3 )?+ (p P <1,0<i <m ,O<j <n ) 


© ¢@:R™ -Dın xR xn -1) _, R™ *" is linear operation, which is defined according 
to the following equation. 
C(p “Pp i J =p : id i -lj i j i „j -1 
Let <m ,1<sj <n 
@ The operator C7 :R™ *7 > R™ -Dxn x Rm x™ -D adjacent to C is calculated 
according to the following formula. 


lT (x )=(p tp >» 


where E R™ *" isthe matrix, whose defining equations are 
8 eq 


p ; P =X i j TX į ay 2! =1,.,m -1,j =1,..,n 
2 _ _ . -_ . = E 
p ; F =X į j X i jari =1,..,,.m ,j =1,..,n 1 
© P œ belongs to an orthogonal projection operator in set C . Therefore, when C = 
B , , »theP 3 can be defined as: 
l „u 
l X ij <l 
P $ _— (X ) j =4X ij l SX ij <u 
u x ij >u 


Based on the notation above, we obtain the dual problem in Eq. (33), and show the relationship 


between the primal and dual optimal solutions. To analyze this relationship, we propose the following 
hypothesis: We assume that ((p 1',p ?)) €P isthe optimal solution to this problem. 
h(p 'p ?°)=-IH ¢ (b -A C(p tp *)I + 
min 


1 


2 
(p ,P 2)EP \|b a À C(p p DI 
The H œ (x )defined as: 


(34) 


H o (x )=x -P ¢ (x ) x €R™™ (35) 


Because objective function h in Eqs. (34) is continuously differentiable and its gradient is defined 


as 
Vip !,p ?)=-2à CT P o (b -à CQ tp *) (36) 
T (x p 4p %2=Tr č (Cp Lp %) x) (37) 
Therefore, problem Eq. (34) can be converted into the following equation. 
2 
min max {Ix -b | 
x EC (p Lp EP a (38) 
T 
+2ATr (cp Lp 3] x )} 


Because the objective function is concave in p !1,p ° and convexinx , the maximum and 


minimum orders can be exchanged. Thus, the following equation is obtained: 


max min {lix —b i 
(p lp EP X EC F (39) 
+2ATr (C(p 1,p 2)” x )} 


This can be rewritten as follows: 
2 
max min {|X -(b -A £ an 2 

o XZ ymin {Ix -( (w 4p a) 


_ _ 1 2 2 2 
lb -à Cp Lp 3È +b Ñ} 


where the optimal solution of the minimum problem is 


x =P ¢ (b =A Cip Lp *)) (41) 
Subsequently, by substituting Eq. (41) into Eq. (40) and omitting the constant term, we obtain the 


|p (40) 


following dual problem: 


P a (> =4 ip Ap J0 =a ip. S A 
max 5 (42) 
Ain =i -A Co Lp UE 
Our goal is to solve the dual problem in Eq. (34), whose gradient is expressed by Eq. (36). 


Therefore, the gradient projection algorithm, which is an effective method for solving the denoising 


problem, is presented. Because the norm of x€ R” ~r is Frobenius norm. For 
(p tip De R™ -Dxrn xy R™ xM -1 the norm is expressed as follows: 
lo Lp D=ylp 1% +ip Ñ (42) 
F F 
L (h)sa ? (43) 
The projection onto set P can be computed simply: For (p lp 2) , projection 


P p (p tp 7) is given by P p (p tp *)=(r },r_ ?). The equations of 


re R ™ *" — are shown below: 
1 
p ij 
max{1,,(p + | )2?+(p ? 3 
1 = ij ij 
> E 1 (42) 

; EF 

max {1, |p : i |+|P ; , p 


2 


re. 
max {1, |p 1 24 (p 2 —_ +92} 
r2? = E E 
I 4 
J p F 
max {1, |p > F |+|P ? ; p 


Substituting the objective function, gradient equation, and maximum Lipschitz constant into Eqs. 
(34) into the gradient projection algorithm introduced above, we obtain the following denoising 
algorithm: 


2.3.4 OS-SART algorithm 


In 1984, Anderson and Kak proposed an improved SART. The SART improves on the ART 
algorithms by calculating the error of all rays passing through a pixel at the same projection angle. The 
ART uses only one ray per iteration, whereas the SART corrects all the rays at the same angle. This is 
similar to lowering the noise introduced by the ART algorithm. Its formula is as follows: 


M-1 
= (k) 
> A Pi ae WinXim 
PicPo k M-I Wy 
ae Wim (30) 
J J 
> Pi EPo Wy 


The OS-SART is a combination of ordered subsets and the SART algorithm, and its formula is as 
follows: 


M-1 
(k) 
pa A Pi E ys WinXim 
ieS, k 
(k+1) = k 


M-1 W; 
Da Wim (31) 
J ay k 
Dns Wi 


where i denotes the i-th ray at a certain projection angle, j denotes the pixel value in the i-th ray, k 


Xx 


denotes the number of iterations, and À is the relaxation factor. 


2.3.4 Proposed algorithm 


Based on the OS-SART iterative algorithm and first-order primal-dual algorithm with total 
variation denoising, we propose a new algorithm for sparse-view NCT 3D image reconstruction (OS- 
SART-PDTYV). The iterative process of this algorithm contains two loops: the outer and inner loops are 
OS-SART for image reconstruction and PDTV for image denoising. 

Table2 OS-SART-PDTV algorithm 


1. While the stop criterion is not met 


2. Stepl: OS-SART reconstruct the image 
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6. Non-negativity constraint; 
If g0 < 0,20 =0 


OS-SART OS-SART 


Step2: PDTV 
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10. While the stop criterion is not met 


11. g“? =(H'H+ (H u+ +g” ~aR"E) 
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16. Output recovery image g+” 


17. end if; 


18. Until the stopping criteria are reached. 


19. Get the final reconstructed image 2o¢_ cuar_ppry 


Il EXPERIMENT 


3.1 Quantitative Evaluation Index 


In this study, we used four image evaluation metrics to quantitatively analyze the mass of the 
reconstructed images for each algorithm. The accuracy of the image reconstruction was quantitatively 
analyzed using the correlation coefficient (CC) metric, which is defined as follows: 


(EO) BY ne) Fine) } 


is the image to be reconstructed, M denotes the voxel number, and @(y) represents 


CES (31) 


where O pue 


the reconstructed value at voxel y . When the reconstructed image was the same as the image to be 


reconstructed, the CC value was 1. 

The ((MSSIM) provides a comprehensive evaluation of the reconstructed image by comparing the 
differences between the image to be reconstructed and the reconstructed image in terms of brightness, 
contrast, and structure. 
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where ¢,7 denote the standard image and the image to be evaluated, respectively; u, and u, are 


the mean values of the image; o, and o, are the standard deviations of the image; ¢,, is the 


covariance of the image; C 1,C 2 andC_ 3 are constants. The mean structural similarity of an 
image is defined as 


L a P 
MSSM = enn Aeon xien (2) 
t=1 
RMSE is defined by the following equations. 
N pa 2 
RMSE = [Zaer (Pn = Pr) (33) 
N 


P,, represents a voxel in the image to be reconstructed, p, represents a voxel in the reconstructed 


sample, and N indicates the number of voxels. A higher RMSE value indicates a larger error between 
the reconstructed images. 

Universal quality image (UQI) is a widely used image evaluation index, which is defined as 
follows: 


UOI = 2cov(u, 4) 244M 


m2 2 m2 2 (33) 
P +P H+H 
where C OV denotes the covariance function, “7, 4 are the means, and p>; p are the 


variances of the reconstructed image and the image to be reconstructed, respectively. The UQI value 
ranges from 0 to 1 and increases with similarity. 


3.2 Digital 3D Shepp—Logan Model Experiment 


In the digital simulation experiment, we compared and analyzed the performance of five 
algorithms: FBP, OS-SART-TV, OS-SART-AwTV, OS-SART-FGPTV, and OS-SART-PDTV by 
reconstructing 3D Shepp—Logan model projection data. The dimensions of the 3D Shepp—Logan model 
were 256x256x256. An iterative algorithm must set reasonable iterative parameters to obtain an ideally 
reconstructed image. The following parameters were used to compare the three iterative algorithms: 


OS-SART-TV à = 2.7, Ni t er gs SaR =80, à py = 200, 
and Ni ter ry =300 . For the OS-SART-AwTV algorithm à =2.5 , 
Niter OS -SART =80 , À awry = 100 , and 
Niter cry =250 . For the OS-SART-FGPTV algorithm, à =2.2 , 
Niter OS . SART = 80 ; MepTry = 0.006 ; 
Niter FOPTV = 150. For the OS-SART-PDTV algorithm, à = 2.2, 
Niter OS -SART = 80 ; AP DTV = 0.0032 : 
Niter PDTV = 300. To verify the denoising capability of the algorithm, the 


following parameters were set according to the noise model of the NCT system: These parameters 


include the photon incident flux 1e 5,m ,;, = 0, and 6 = 10. 

Fig.2 shows the relationship between the RMSE of the reconstructed images using the OS-SART- 
PDTV algorithm and the number of iterations for different sparse views (28 and 38 projection views), 
showing that the OS-SART-PDTV algorithm can converge to a stable solution after a certain number of 


iterations. 
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Fig.2 RMSE versus iteration steps for the OS-SART-PDTV algorithm from different sparse views: (a) 
28; (b) 38. 


Fig. 3 shows the reconstruction of images with different numbers of projected views using the five 
algorithms. As illustrated in Fig. 3, the smoothness and sharpness of the images reconstructed using the 
five algorithms improved rapidly as the number of projection views increased. Compared to the other 
three iterative algorithms, the reconstructed images of the FBP algorithm showed severe bar artifacts. 
The reconstructed images from the four iterative algorithms are smoother and have sharper boundaries 
than those from the FBP algorithm. As shown in the enlarged ROI in Fig. 3, the OS-SART-TV, OS- 
SART-AwTV, and OS-SART-FGPTV algorithms can effectively reduce artifacts; however, part of the 
edge structure is not exactly reconstructed. In contrast, the novel algorithm can not only maintain good 
reconstruction performance but also reconstruct finer structures. 
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Fig.3. Sliced images at z=130 for four algorithms reconstructing reconstructed images with different 
numbers of projection views. 

In addition, we compared and analyzed the error images of the reconstructed images using 
different algorithms that reflected the difference between the pixel values of the reconstructed and 
reference images. According to the error image in Fig. 4, the image reconstructed using the FBP 
algorithm lost many features and produced many artifacts, indicating that the algorithm is not suitable 
for sparse-view NCT 3D reconstruction. The other four iterative algorithms performed better at 
suppressing artifacts and reconstructing detailed structural information. The OS-SART-PDTV 
algorithm has much less detail loss in the error image than the other algorithms. 

We plotted the profile of the reconstructed images for each algorithm in Fig. 5 to further evaluate 
the performance of each algorithm. The number of projection views from left to right was 28, 38, and 
38 +, respectively. Here, 38+ indicates that the noise model was added to the 3D Shepp-Logan model. 
As shown in Fig. 5, there was a large fluctuation in the FBP profile line with the largest deviation from 
the true pixel value. Although the OS-SART-TV, OS-SART-AwTV, and OS-SART-FGPTV algorithms 
performed better than the FBP algorithm, there was still a gap between their amplitudes and true pixel 
values. The profile of OS-SART-PDTV matches the reference value best and is closest to the true value. 
This also indicates that the OS-SART-PDTV algorithm is well-suited for sparse-view NCT 3D 
reconstruction. 
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Fig.4. Slice of error images at z=130 
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Fig.5. Horizontal profiles of the 3D Shepp—Logan model at z =130" slices. (a) 28; (b) 38; (c) 38+. 
As shown in Fig. 6. We used four image evaluation metrics to quantitatively analyze the 
reconstructed images from each algorithm. First, the evaluation indices of the images reconstructed 
using noiseless projection data were analyzed. From the CC, MSSIM, and UQI histograms shown in 
Figure 6, the values of the three metrics for the reconstructed images from each reconstruction 
algorithm gradually increased as the number of projection views increased. These three metrics exhibit 
similar characteristics. 
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According to the RMSE histogram in Fig. 6c, the reconstructed images of OS-SART-PDTV had 
the smallest RMSE values compared to the other algorithms for the same projection views. The RMSE 
of the reconstructed image for each algorithm gradually decreased as the number of projection views 
increased. The RMSE value of the images reconstructed using the OS-SART-PDTV algorithm at 28 
projection views was 0.03352, indicating that the reconstructed image was closest to the real image and 
that the quality of the reconstructed images was higher. As shown in the UQI histogram in Fig. 6d, the 
UQI values of the reconstructed images for each algorithm gradually increased as the number of 
projected views increased. The trend of the UQI was the opposite to that of the RMSE. The 
reconstructed images from the OS-SART-PDTV algorithm had the highest UQI values for the same 
number of projection views, indicating that it outperformed the other algorithms. The UQI value of the 
image reconstructed using the OS-SART-PDTV algorithm was 0.99039 for the 28 projection views. We 
then quantitatively analyzed four evaluation metrics for the images reconstructed using data containing 
noisy projections. From the CC, MSSIM, and UQI histograms, it can be seen that the OS-SART-PDTV 
algorithm has the largest evaluation metrics for these three reconstructed images when the same 
number of noise-containing projections are reconstructed. According to the RMSE histogram, the 
reconstructed image from the OS-SART-PDTV algorithm had the lowest RMSE value. Therefore, 
based on the quantitative and visual analyses of each algorithm, the OS-SART-PDTV algorithm has the 
highest-quality reconstructed images with the same number of projection views. 
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Fig.6. Histogram of reconstructed image evaluation metrics for the 3D Shepp-Logan model. (a) CC; (b) 
MSSIM; (c) RMSE; (d) UQI. 


3.3 Digital Head Model Experiment 
We further analyzed the performance of these five algorithms by reconstructing a digital head 


model. The reconstructed images from the five algorithms used for reconstructing different projection 
view numbers are shown in Fig.7. Based on the reconstructed images, the smoothness and sharpness of 


the reconstructed images of all five algorithms improve significantly as the number of projected views 
increases. Compared to the other four algorithms, the reconstructed images of the FBP algorithm show 
severe bar artifacts. The reconstructed images from these four iterative algorithms were smoother and 
had cleaner image boundaries than those from the FBP algorithm. As shown in the enlarged ROI in 
Fig.7, the OS-SART-TV, OS-SART-AwTV, and OS-SART-FGPTV algorithms are effective in reducing 
artifacts, but the edge structure of the image is not accurately reconstructed. In contrast, the proposed 
OS-SART-PDTV algorithm can not only reduce artifacts, but also reconstruct more detailed structural 
information. 

To further evaluate the performance of the algorithm, the contours of the reconstructed images for 

each algorithm are plotted in Fig. 8. As shown in Fig.8, the numbers of projected views from left to 
right are 28, 38, and 38+, where 38+ indicates that noise has been added to the projection data. As 
shown in Fig.8, there was a large fluctuation in the FBP profile, which had the largest gap from the true 
pixel values. Although the OS-SART-TV, OS-SART-AwTV, and OS-SART-FGPTV algorithms 
performed better than the FBP algorithm, there was still a large deviation between the actual and true 
pixel values. The profile of OS-SART-PDTV matches the reference value best and is closest to the true 
value. This also indicates that the OS-SART-PDTV algorithm is well-suited for sparse-view NCT 3D 
reconstruction. 
According to the RMSE histogram of the reconstructed images from each algorithm, the OS-SART- 
PDTV algorithm had the smallest RMSE of the reconstructed images. As the number of projected 
views gradually increased, the RMSE of the reconstructed image for each algorithm gradually 
decreased. When the number of projected views is 28, the RMSE value of the OS-SART-PDTV 
algorithm-reconstructed image is 0.02404, which indicates that the reconstructed image of this 
algorithm is the closest to the real image, and the quality of the reconstructed image is higher. 

The UQI histogram in Fig. 9d shows that the UQI values of the reconstructed images of each 
algorithm increased as the number of projected views increased. The UQI trend was opposite to that of 
the RMSE. The OS-SART-PDTV algorithm reconstructs the image with the largest UQI value for the 
same number of projected views, indicating that this algorithm outperforms the other three iterative 
algorithms. When the number of projected views was 28, the UQI of the image reconstructed by PDTV 
algorithm was 0.98594. We then quantitatively evaluated images containing noisy projections 
reconstructed using the four algorithms. The OS-SART-PDTV algorithm reconstructs the image with 
the largest CC, MSSIM, and UQI values for the same number of projected views. According to the 
RMSE histogram, the RMSE value of the reconstructed image obtained using OS-SART-PDTV was 
the smallest. Therefore, by quantitative and visual analyses of each algorithm, OS-SART-PDTV 
exhibited superior performance in sparse view reconstruction. 

To analyze the relationship between the OS-SART-PDTV algorithm and regularization parameters, 
we present the reconstructed images with different regularization parameters in Fig.10. Based on the 
reconstructed images, an optimal image was obtained by adjusting the regularization parameters when 
reconstructing the head model. According to the reconstructed images with different regularization 
parameters, it can be seen that the reconstructed image quality reaches the best when the regularization 
parameters reach a certain value as the reconstruction parameters gradually increase. As the 
regularization parameter continued to increase, the quality of the reconstructed image decreased. 
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Fig.7. Sliced images at z=70 for four algorithms reconstructing reconstructed images with different 
number of projection views. 
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Fig.8. Horizontal profiles of the digital head model at z=250" slices. (a) 28; (b) 38; (c) 38+. 
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Fig.9. Histogram of reconstructed image evaluation metrics for the digital head model. (a) CC; (b) 
MSSIM; (c)RMSE; (d)UQI. 
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Fig.10. Reconstructed images with different regularization parameters for the OS-SART-PDTV 
algorithm 38+ 


3.3 Real neutron projection experiment 


To further validate the performance of this novel algorithm in practical applications, the 
performance of the algorithm was verified using projection data from an NCT-based clock model 
(Fig.11) provided by Burkhard Schillinger of the Technical University of Munich, Germany [37]. 
Schillinger et al. provided 201 neutron projection images captured at equal angles in the range of 0 °to 
180°. They also provided two dark field and two open field background images. The background 
images were obtained using the CCD camera with the neutron beam turned off, which contained 
readout noise as well as electron noise caused by dark currents. Open-field images were obtained using 
a CDD camera with neutrons turned on and without the samples. 
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Fig.11. Clock Model 

As shown in Fig.12, the artifacts of the reconstructed images by the FBP algorithm gradually 
increase as the number of projected views decreases. When the number of projected views is 45, the 
four iterative algorithms outperform the FBP algorithm in suppressing artifacts and reconstructing fine 
structure information. Although the OS-SART-TV, reconstructs images with fewer artifacts, the edges 
of the reconstructed images are too smooth, resulting in the loss of detailed structural information of 
the reconstruction images. Despite OS-SART-AwTV and OS-SART-FGPTV algorithms can also 
reconstruct a lot of image details, their reconstruction results are poor compared with the OS-SART- 
PDTV algorithm that reconstructs the same number of sparse views. Visually, the OS-SART-PDTV 
algorithm outperforms other algorithms in terms of noise removal, artifact suppression, and retention of 
detailed structural information. 

The performances of the three iterative algorithms were quantitatively analyzed using four image- 
evaluation metrics, as shown in Fig.13. In Fig.13, the RMSE of the image reconstructed by the OS- 
SART-PDTV algorithm is the smallest, indicating the minimal error between the reference image and 
the reconstructed image, whereas the image reconstructed by the OS-SART-PDTV algorithm has the 
highest CC, MSSIM, and UQI values, indicating that it is most similar to the reference image. In 
summary, the OS-SART-PDTV algorithm performed well in denoising, suppressing artifacts, and 
reconstructing fine structures when reconstructing the projection data of sparse NCT scans, indicating 
that this novel algorithm can be applied to the 3D reconstruction of sparse-view NCT. 
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Fig.12. Sliced images at z=201 for four algorithms reconstructing reconstructed images with different 


number of projection views. 
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Fig.13. Histogram of reconstructed image evaluation metrics for the clock model. (a) CC; (b) MSSIM; 
(c)RMSE; (d)UQI. 
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Fig.14. The profile of reconstructed images by FBP, OS-SART-TV, OS-SART-AwTV, OS-SART- 

FGPTV, and OS-SART-PDTV at z=201. The number of projection views is (a) 45; (b) 90; (c) 135. 

Fig.14 presents the profiles of the reconstructed images. According to the profile, the OS-SART- 
PDTV algorithm has an advantage over the other algorithms, and its profile is the closest to the 


reference. 
IV DISCUSSION 


Herein, we proposed a 3D reconstruction algorithm for sparse-view NCT. Because the analytical 
algorithm (FBP) has inherent defects that make it unsuitable for sparse-view reconstruction, we used an 
iterative algorithm applied to the 3D reconstruction of sparse-view NCT. Furthermore, the iterative 
algorithm can embed a priori information to regularize the image. Therefore, the TV-based iterative 
algorithm is widely used in sparse-view CT image reconstruction because it is superior to other 
algorithms in artifact suppression and denoising. By reconstructing the projection data of the 3D 
Shepp-Logan, head, and clock models, the OS-SART-PDTV algorithm outperformed the other 
algorithms in both quantitative and visual analyses. Therefore, OS-SART-PDTV is suitable for sparse- 


view NCT. In addition, deep learning, an important research direction in the field of image 
reconstruction, has led to many significant research results in sparse-view reconstruction, such as in 
[38], using short C-arm CT. scans. This study proposed the DRONE model, which uses a codec 
network in the embedding module to extract depth features in the data and image domains, combined 
with the Wasserstein distance generation adversarial network to maintain details and features in the 
image domain. It combines data residual and image residual networks in the refinement module to 
restore fine structure features from the output of the embedding module. According to the CS iterative 
reconstruction model, the depth prior of the data and image domains is regularized, thus ensuring the 
robustness of the DRONE network and making the DRONE model good at retaining image edge 
information, recovering image features, and reconstructing accurate image information in practical 
applications. 

Generally, in practical NCT scanning systems, the projected image is affected by photon statistical 
noise and electrical noise. Therefore, several iterative reconstruction algorithms cannot achieve 
satisfactory results in NCT scanning systems. To further test the performance of our algorithm, we 
reconstructed the projection data obtained from the NCT scan-clock model using the OS-SART-PDTV 
algorithm. Based on the results of the reconstructed neutron projection data, OS-SART-PDTV 
exhibited a good performance in denoising, suppressing artifacts, and reconstructing fine structures. 
According to the results of the reconstruction of the 3D Shepp-Logan model, the OS-SART-PDTV 
algorithm converged monotonically to a stable solution (Fig. 3). 

Iterative reconstruction algorithms generally need to optimize several parameters to obtain a stable 
solution. This is a common problem encountered by all iterative reconstruction algorithms. The OS- 


SART-PDTV algorithm must optimize the following parameters: A os -SART ; 
N os -SART „À pPprv and N ppry . We found that the quality of 
the reconstructed images was significantly improved when À os -SART varied within a 


small range. In other words, the OS-SART-PDTV algorithm was more sensitive to this parameter. 
N os -SART is usually set in the range of 50—100 because the improvement rate of the 
image quality gradually decreases as the number of iterations increases, and after reaching a certain 
number of iterations, the image quality no longer improves. According to our experimental results, the 
OS-SART-PDTV algorithm obtained a better reconstructed image after 50 iterations. In addition, the 
two parameters À ppry and N ppry must be tuned several times to obtain high- 
quality reconstructed images. When reconstructing the projection data of a new sample, the two 


parameters À ppry and N ppry must be readjusted. 


V CONCLUSION 


NCT has proven to be a highly effective noninvasive measurement method and has been widely 
applied in nuclear engineering, thermodynamics, and other fields. However, because of the complexity 
of the NCT application environment, obtaining complete projection data is challenging. The selection 
of an efficient sparse-view reconstruction algorithm is directly related to the system imaging quality 
and reliability of the measurement results. To this end, we propose the OS-SART-PDTV algorithm for 
sparse-view NCT 3D. The OS-SART-PDTV algorithm is divided into two steps: OS-SART for image 
reconstruction, and PDTV for solving the total variance using a first-order primal-dual algorithm. 
Based on the reconstruction results of the projection data obtained from the NCT scan clock model, the 
OS-SART-PDTV algorithm has a significant advantage over the other algorithms in preserving edge 
structure, denoising, and suppressing artifacts. Therefore, the OS-SART-PDTV algorithm has great 
potential for application in NCT 3D reconstruction. 
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